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Dose-finding studies are frequently conducted to evaluate the ef- 
fect of different doses or concentration levels of a compound on a re- 
sponse of interest. Applications include the investigation of a new 
medicinal drug, a herbicide or fertilizer, a molecular entity, an envi- 
Ph ' ronmental toxin, or an industrial chemical. In pharmaceutical drug 

•^r ' development, dose-finding studies are of critical importance because 

of regulatory requirements that marketed doses are safe and pro- 
C^ ' vide clinically relevant efficacy. Motivated by a dose-finding study in 

^ moderate persistent asthma, we propose response-adaptive designs 

addressing two major challenges in dose-finding studies: uncertainty 
about the dose-response models and large variability in parameter 
estimates. To allocate new cohorts of patients in an ongoing study, 
we use optimal designs that are robust under model uncertainty. In 
QQ ' addition, we use a Bayesian shrinkage approach to stabilize the pa- 

QQ , rameter estimates over the successive interim analyses used in the 

l/^ ' adaptations. This approach allows us to calculate updated parameter 

l»^ ^ estimates and model probabilities that can then be used to calculate 

f— V ■ the optimal design for subsequent cohorts. The resulting designs are 

hence robust with respect to model misspecification and addition- 
ally can efficiently adapt to the information accrued in an ongoing 
study. We focus on adaptive designs for estimating the minimum 
effective dose, although alternative optimality criteria or mixtures 
^^ ■ thereof could be used, enabling the design to address multiple objec- 

$_^ ' tives. In an extensive simulation study, we investigate the operating 

characteristics of the proposed methods under a variety of scenarios 
discussed by the clinical team to design the aforementioned clinical 
study. 
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1. Introduction. Dose-finding studies liave several challenges in com- 
mon. First, they usually address two distinct objectives, which lead to differ- 
ent requirements on the study design [Ruberg (1995), Bretz et al. (2008)]: 
(i) assessing evidence of a drug effect, and (ii) estimating relevant target 
doses. Second, the form of the dose-response relationship is unknown prior 
to the study, leading to model uncertainty. This problem is often underesti- 
mated, although ignoring model uncertainty can lead to highly undesirable 
effects [Chatfield (1995), Draper (1995), Hjorth (1994)]. Third, data from 
dose-finding studies are usually highly variable. This issue is of particular 
importance in pharmaceutical drug development, because sample sizes are 
kept to a minimum for ethical and financial reasons. It is therefore critical to 
develop efficient dose-finding study designs that use the limited information 
as efficiently as possible, while addressing the above challenges. 

Many approaches have been proposed in the optimal design literature to 
distribute patients efficiently with regard to given study objectives; see Wu 
(1988), Fedorov and Leonov (2001) and King and Wong (2004), among many 
others. However, most of this work has concentrated on an assumed fixed 
dose-response model. As there is typically considerable model uncertainty 
at the planning stage of a dose-response study, these methods have limited 
practical use. Based on concepts introduced by Lauter (1974) [see also Cook 
and Wong (1994), Zhu and Wong (2000, 2001), Biedermann, Dette and Pe- 
pelyshev (2006)] , Dette et al. (2008) investigated model-robust designs that 
provide efficient target dose estimates for a set of candidate dose-response 
models, rather than for a single dose-response model. However, their de- 
signs require knowledge about the unknown parameters associated with the 
anticipated dose-response models as well as the prespecification of model 
probabilities. 

A natural remedy is to investigate response-adaptive designs (adaptive 
designs, in short) with several cohorts of subjects. After each stage the ac- 
cumulated information of the ongoing study is used to update the initial 
information of the underlying model parameters and model probabilities, 
which in turn is used to calculate the design for the subsequent stage(s). 
Several adaptive designs have been developed for this problem; see, for ex- 
ample, Miller, Guilbaud and Dette (2007) and Dragalin, Hsuan and Pad- 
manabhan (2007) for recent approaches using optimal design theory, or Zhou 
et al. (2003), Miiller et al. (2006), and Wathen and Thall (2008) for recent 
Bayesian adaptive designs. Dragalin et al. (2010) performed an extensive 
simulation study that compared five different adaptive dose-finding meth- 
ods. 

In this paper we propose adaptive designs addressing the three major 
challenges described above: multiple study objectives, model uncertainty 
and large variability in the data. For this purpose we use the model-robust 
designs proposed by Dette et al. (2008) together with a Bayesian shrinkage 
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approach to stabilize the parameter estimates, especially in the early part of 
a study. This allows one to calculate parameter estimates as well as model 
probabilities that can then be used to calculate model-robust designs for 
the subsequent stage(s) of the study. The resulting designs are robust with 
respect to model misspecification and additionally adapt to the continuously 
accrued information in an ongoing study. We focus on adaptive designs for 
estimating the minimum effective dose (MED), that is, the smallest dose 
achieving a clinically relevant benefit over the placebo response. However, 
alternative optimality criteria or mixtures of optimality criteria could be 
used, enabling the design to address multiple objectives. 

2. Asthma dose-finding study. The research for this article was moti- 
vated by a Phase II dose-finding study for the development of a new pharma- 
ceutical compound in asthma. This was a multi-center, randomized, double- 
blind, placebo controlled, parallel group study in patients with moderate 
persistent asthma, who were randomized to one of seven active dose levels 
or placebo. The primary endpoint was change from baseline in a lung func- 
tion parameter (forced expiratory volume in 1 second, FEVi) after 28 days 
of administration, scaled such that larger values indicated a better outcome. 
The objective of the trial was to evaluate the dose effects over placebo for 
the primary endpoint and to assess whether there was any evidence of a drug 
effect. Once such a dose-response signal had been detected, one would sub- 
sequently estimate relevant target doses, where the primary focus was on 
estimating the MED. 

Based on discussions with the clinical team, a homoscedastic normal 
model was assumed for the primary endpoint with a standard deviation 
of 350 ml, a placebo effect of 100 ml and a maximum treatment effect of 
300 ml within the dose range [0, 50] under investigation. The available doses 
were (= placebo), 0.5, 1.0, 2.5, 5, 10, 20 and 50. The clinically relevant 
benefit over the placebo effect was set to 200 ml. That is, an increase in 
treatment effect of less than 200 ml over the observed placebo response was 
considered to be clinically irrelevant. Furthermore, all dose levels within the 
investigated dose range were considered safe based on previous studies, so 
that efficacy was of primary interest. 

Because this study was conducted early in the drug development pro- 
gram, limited information about the dose-response shape was available at 
the planning stage. A set of candidate dose-response models was derived 
before starting the study; see Table 1 and Figure 1 for the full model spec- 
ifications (including a preliminary specification of the model parameters). 
An increase of the dose-response curve in the lower part of the investigated 
dose range was considered likely, so two concave increasing models (Emaxi , 
Emax2) were included in the model set. In addition, 5-shaped (Logistici), 
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unimodal (Beta) and convex (Logistic2) models were included in the can- 
didate model set to robustify the statistical analysis with respect to model 
uncertainty. We refer to Pinheiro, Bornkamp and Bretz (2006) for details on 
the use of candidate models in dose-response studies and the elicitation of 
best guesses for the model parameters. 



Model 



Table 1 

Candidate dose response models as a function of dose d, where 

B{a,b) = {a + by+''/{a''b'') 



Full model specification 



Model parameters True MED 



Beta 00 + 6(iS(6'2, 6'a)(rf/60)®2 (1 - d/60)* 
Emaxi eo + eid/{e2 + d) 

Emax2 eo + eid/{e2 + d) 

Logistici 00 + 0i/{l + exp[(6'2 - d)/e3]} 

Logistica Oo + 0i/{l + exp[(6»2 - d)/es]} 
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Fig. 1. Graphical display of the dose-response models in Table 1. Open dots denote the 
potential responses at the seven active dose levels and placebo available in the study. Dotted 
horizontal lines indicate the clinical relevance threshold on top of placebo response and 
dotted vertical Imes the resulting MED. 
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Given the information and constraints above, the chnical team was faced 
at the planning stage with several remaining key questions on the study 
design: 

(A) Should an adaptive design be employed at all or would a nonadaptive 
design be sufficient? 

(B) If the decision was to employ an adaptive design, how many interim 
analyses should be conducted? 

(C) How many dose levels should be included in the study, that is, are 
all seven active dose levels from above needed? 

(D) If not all active dose levels were needed, which of them should then 
be investigated? 

In addition to these statistical questions, many further considerations were 
discussed by the clinical team: adaptive designs require more logistical effort 
to set up the repeated data collection and cleaning/ analysis processes than 
nonadaptive designs; including all seven active doses in the study would 
pose serious challenges to the drug manufacturing and supply departments, 
especially if the allocation changed during the study course; and how to 
ensure trial integrity and validity. In the following we focus on the statistical 
questions and describe the proposed methodology for the study. 

3. Methodology. Assume k distinct dose levels di, . . . ,dk, where di = 
denotes placebo. Let n^ patients be allocated to dose di and -/V = X]j=i"'i- 
The vector of allocation weights is denoted by w = {wi, . . . , Wk)', where Wi = 
Ui/N. Let further Yij r^j ]S[{f(di,0),a'^) denote the observation of patient 
J = 1, . . . , rij at dose di,i = 1, . . . ,k, where the dose-response model /(•) is 
parameterized through the parameter vector 6 and N{fi,a'^) denotes the 
normal distribution with mean fi and variance o"^. 

Most dose-response models used in practice, including those in Table 1, 
can be decomposed as 

(3.1) fid,e) = eo + 6if{d,e''), 

where 9 = (6*0, 6*1, 61°')' = (e^o, . . . ,6*^)'. The parameters 6* = (6*0, 6*1)' enter 
the model function / linearly and determine its location and scale, while /" 
is typically a nonlinear function that determines the shape of the model 
function / through the parameters 9 . 

The minimum effective dose producing a clinically relevant effect A over 
the placebo response is defined as 

(3.2) MED= min {f{d,9)>f{di,9) + A}, 

where we assume that a beneficial effect is associated with larger values 
of the response variable. Note that the MED may not exist, as no dose in 
((ii,(ifc] may produce an improvement of A compared with placebo. 
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3.1. Robust designs for MED estimation. Given a function /^, it follows 
from (3.2) that the MED (provided it exists) is a solution to 

(3.3) ^0 + ^i/°(0, 0°) + A = 00 + OifiMEB, 0°). 

Consequently, MED = 6(0) = h^{f{0,e^) + A/9i), where h^{x) = mi{z\ 
f^{z) > x} denotes the (generalized) inverse of the function f^ with respect 
to the variable d. Standard asymptotic theory for nonlinear models [Se- 
ber and Wild (1989)] yields that the maximum likelihood (ML) estimate 6 
is approximately multivariate normal distributed with mean vector 6 and 



,2 



covariance matrix ^M ^(0,w), where M.{6,-w) = '^-^^Wig{di,9)g{di,d)' 

denotes the information matrix and g{d, 6) = »g' ' the gradient of the 
dose-response model / with respect to 0. It follows from the (5-method 
[see van der Vaart (1998)] that the MED estimator based on 0, MED = 
b{0), is asymptotically normally distributed with mean h{0) and variance 
y(0,w) = ^V6(0)'M^i(0, w)V6(6>), where V5(0) = ^. Hence, minimi- 
zing V{0, w) with respect to w G S'^ = {w| X]i=i ""^i = 1; w > 0} results in op- 
timal designs that minimize the approximate variance of MED . This design 
criterion has also an appealing decision theoretic justification: The asymp- 
totic normal distribution of MED approximates the posterior distribution 
of the MED in a Bayesian model framework. Hence, minimizing the log- 
variance of MED is equivalent to minimizing the (approximate) Shannon 
entropy of the posterior distribution of the MED [Chaloner and Verdinelli 
(1995)]. 

In principle, the above optimization could be done with respect to the 
number and choice of doses and their corresponding allocation ratios [Dette 
et al. (2008)], but, in practice, manufacturing constraints often determine the 
available doses, as it was the case in the asthma study from Section 2. In the 
following we thus restrict the optimization to the weights w for prespecified 
doses di, . . . ,dk- 

The true dose-response function / is unknown and optimal designs are ty- 
pically not robust with respect to model misspecification [Dette et al. (2008)]. 
In the following we assume a set of M candidate models fm{d, 9m) = Oom + 
dim f mid: ^m); rn = 1,. . . ,M, such as those described in Table 1. We "inte- 
grate" the design criterion conditional on model m with respect to the model 
probabilities Om- Hence, using the design criterion J2m=i ^m log(Kra(^m, w)) 
or, equivalently, 

M 

(3.4) ^M=ll{Vm{em,^)r"' 

m=l 

leads to designs that are robust with respect to model misspecification, where 
V^(0m)W) denotes the variance of the estimate for the MED in the inih. 
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model {m = 1, . . . , M). Note that because of taking logarithms above, there is 
no need to standardize the individual model variances. Otherwise this would 
be necessary to avoid that some models dominate the design criterion [the 
Vm{Om,'vv) can be quite model dependent and differing in size]. However, the 
numerical calculation of robust designs using the criterion (3.4) requires the 
knowledge of 0^ and Om, m = 1, . . . , M. In the following sections we describe 
how the initial best parameter guesses can be updated during an ongoing 
study such that subsequent stages can be redesigned based on the updated 
estimates for 6m and am,m = 1, . . . ,M; see Section 3.3 for a description of 
the complete procedure in an algorithmic form. 

3.2. Updating of model parameters and weights. Reliably estimating the 
parameters 6i, . . . , 9m is a challenging problem, particularly in early stages 
of a study. ML estimates for these parameters are typically highly variable, 
and may even not exist without imposing bounds on the parameter space. 
One way of stabilizing estimates is to use a shrinkage approach based on, 
for example, penalized maximum likelihood or maximum a-posteriori (MAP) 
estimates. Here, one optimizes the log-likelihood function plus a term which 
determines the prior plausibility of the parameters (the log prior distribu- 
tion). The estimate is then a compromise between the information con- 
tained in data and the prior distribution. This stabilizes the estimates in 
early stages due to the shrinkage toward a priori reasonable values. In later 
stages the shrinkage effect decreases because the log prior remains constant 
while the log likelihood receives more weight with increasing sample sizes. If 
a completely flat prior distribution is used, standard ML and MAP estima- 
tion coincide, so that using nonuniform priors is desirable. We discuss the 
choice of nonuniform priors in more detail further below. 

Apart from stable parameter estimates 6i, ... ,6m for the dose-response 
models, one needs to update the model probabilities ai, . . . , om at an interim 
analysis. We propose using a probability distribution over the different dose- 
response models and evaluating the posterior probabilities for each model 
after having observed the data; see, for example, Kass and Raftery (1995) 
for a detailed description of posterior probabilities and Bayes factors. These 
posterior model probabilities can then be used in the design criterion (3.4). 
A computationally efficient approach to approximate the posterior model 
probabilities is the Bayesian information criterion (BIG). However, previous 
simulation studies in the context of dose-finding studies showed that the BIG 
approximation frequently favors too simplistic models for realistic variances 
and sample sizes [see Bornkamp (2006)]. Other approximate methods, such 
as fractional Bayes factors, or intrinsic Bayes factors [see O'Hagan (1995) or 
Berger and Pericchi (1996)], either depend on arbitrary tuning parameter 
values or are computationally prohibitive. Thus, for each model we will 
use the exact posterior probabilities resulting from the prior distributions 
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assumed for the MAP estimation. In our case, these probabihties can be 
calculated using efficient numerical quadrature without the need to resort 
to computationally expensive Markov chain Monte Carlo techniques. In the 
remainder of this section we provide details on the prior elicitation and the 
calculation of posterior probabilities. 

3.2.1. Selection of prior distributions for 6^. We utilize the factoriza- 
tion in (3.1) to derive a prior distribution for (^omj^imj^rru'^^)- If ^m were 
known, the nonlinear models would reduce to a linear model. It is therefore 
reasonable to use for a given 0^ the conditionally conjugate normal-inverse 
gamma (NIG) distribution 

p{e*^, a'\el) cc {a')-(^+^y^ eM-{iO*m " t^Y^-HO*^ - /x) + a}/i2a')] 

for (0^,0"^) [see O'Hagan and Forster (2004)], where the parameter 6* is 
defined after equation (3.1), a,u >0, /i € M^ and V G R^^^ denotes a pos- 
itive definite matrix. The NIG distribution marginally induces a bivariate 
t-distribution for 6^ with z/ degrees of freedom, finite mean /x and covari- 
ance matrix a/(i^ — 2)V, provided that i/ > 2. The marginal prior distribution 
for (T^ is given by an IG distribution with mode a/(i^ + 2), mean a/(i^ — 2) 
and variance 2a^/{(z^ — 2)^(z^ — 4)}. It has a finite mean when u > 2 and 
a finite variance when ly > 4. 

To set up the NIG distribution for (0J^,(T^), one can employ available 
information about the placebo effect, the maximum treatment effect and 
the standard deviation. For example, one can choose the marginal bivariate 
t-distribution for 0'^ (conditional on 6^) such that the desired mean and 
covariance are achieved for the placebo effect and the maximum effect of the 
underlying dose-response model. When the linear parameters 0J^ cannot be 
interpreted as placebo and maximum effect, one can use a suitable trans- 
formation to achieve the desired moments. Then, one can adjust a and d so 
that the marginal distribution of o"^ achieves the desired mode. An attrac- 
tive choice is to use v = A, leading to a prior with infinite variance for o"^ 
and a heavy tailed marginal prior for 6^. 

For the nonlinear parameters 0^, we propose selecting suitable bounds 
for the parameters and then eliciting a bounded prior distribution. This is 
typically not difficult, as the interpretation of the nonlinear parameters is 
straightforward, and excessively large parameter values usually correspond 
to a priori unlikely model shapes. We propose using a scaled beta distribu- 
tion B{a,P) with mode equal to the initial parameter guesses. The curvature 
of the prior determines the amount of shrinkage that one is willing to employ 
for the MAP estimates. In the simulations we used the sum S = a + f3 as 
a measure of curvature with S" > 2 to ensure unimodality of the beta distribu- 
tion. Note that already relatively small values of S, such as /S = 10 or 5 = 20, 
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Fig. 2. Beta priors on [0.5,75] with mode 20 and different S values. 

lead to strong shrinkage effects; see Figure 2 for an illustration of the 62 pa- 
rameter in the Emaxi model, where the initial parameter guess is 20. For 
dose-response models with more than one nonlinear parameter, we repeat 
this procedure for all parameters and assume independence among them. 

For selecting prior model probabilities, it is convenient to use a uniform 
distribution across the models unless some models are deemed a priori more 
plausible than others. 



3.2.2. Calculation of posterior probabilities. Let y denote the data avail- 
able at an interim analysis and p{y\Om,m) the likelihood under model m 
with corresponding prior distribution p(9m\'m) and prior model probabili- 
ty p[m). Then the marginal likelihood is given by 

p{m\y) (xp{m) / p{y\em,cr'^ ,m)p{Om,cr^\m) d{e,n,cr'^) 

(xp{m) / p{y\em,cr'^,rn)p{e*^,a'^\e^,m)d{e*^,a^)p{e'l^\m)del^. 

The inner integral in the last equation is the product of a likelihood and 
a conjugate prior distribution. One can hence reduce the integration to 



(3.5) 



p{m\y) (xp{m) / p{y\e^^,m)p{6%^\m)de^,^^ 
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where p{y\6^,m) now denotes the integrated hkehhood. In our apphca- 
tions, the integral (3.5) is one- or two-dimensional over a bounded region 
and hence straightforward to calculate numerically. This allows us to calcu- 
late the marginal likelihoods efficiently, without resorting to Markov chain 
Monte Carlo calculations; see Section 3.4 for details. The posterior model 
probabilities p{m\y) can be obtained by normalizing the marginal likelihoods 
(multiplied by the prior model probabilities). 

We use the maximum 6^ of the marginal posterior p{y\0^^,m)p{6^\m) 
as an estimate of 6^. Conditional on this value, we use the maximum 6^ of 
p{6*^\6 ^,y ^m) as an estimate for 0*^. Therefore, the overall estimate of the 

parameter 9^ is given by 0^ = i^rm^m)- This is a slight variation of the 
MAP approach described above, but reduces further the computational ef- 
fort, as it reuses the calculations from the integration to obtain the marginal 
likelihoods. 

3.3. Main algorithm. We now summarize the complete response-adaptive 
dose-finding design in algorithmic form. 
Before trial start: 

(1) Select a starting design using either a balanced allocation across the 
available doses or an unbalanced allocation based on optimal design consid- 
erations. 

(2) Select candidate dose-response models fm{d,9m)- 

(3) Conditional on 6^, calculate a NIG prior distribution for (0^, '^^) 
based on "best guesses" for the placebo effect, the maximum treatment effect 
and o"^ (together with suitable variability assumptions for both parameters). 

(4) Choose "best guesses" for the nonlinear parameters 6^ and select the 
parameter S. 

(5) Choose prior model probabilities p{m) for the different dose-response 
functions. 

At interim analysis: 

(1) Calculate posterior model probabilities 

(3.6) p{m\y) oc p{m) p{y\6m,m)p{6m.\'m)dejn- 

Exploiting the conjugacy properties of the NIG distribution, this reduces to 
one- or two-dimensional integrals; see Section 3.4 for computational details. 

(2) For each dose-response model, estimate 6„i by using the maximum of 
p{y\6^^,'m)p{6^^\m), where the abscissas calculated in step 1 can be reused. 

Conditional on this value, use the maximum ofp{6'^\6^,y, m) as an estimate 
for 0^ to obtain 6^ 
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(3) Plug the obtained parameter estimates 9m into (3.4) and set am = 
p{m\y). Then, minimize with respect to w G 5, where S = {w € S'^jw = 
(noid + iVncxtWnoxt)(A'oid + ^^nGxt)"^ Wncxt G S''}- Here, Hold denotes the vec- 
tor of sample sizes per dose and A'^Qid the total sample size until the current 
interim analysis. Further, A^noxt denotes the sample size and Wnoxt £ § the 
design weights for the next cohort of patients. We therefore optimize the 
design for the next stage taking into account the patient allocation until the 
current interim analysis; see Section 3.4 for computational details. 

(4) Allocate the next cohort of patients according to ■Wncxt by applying 
an appropriate rounding technique, such as described in Pukelsheim (1993), 
Chapter 12. 

Note that the Bayesian approach is used here for design adaptation pur- 
poses. The final analysis may or may not be done using a fully Bayesian 
approach. The development of the Bayesian design methodology above is 
motivated by the MCP-Mod methodology described in Bretz, Pinheiro and 
Branson (2005) to address model uncertainty. This method requires prior 
estimates for the placebo effect, the maximum treatment effect, and a'^ at 
the design stage and "best guesses" of the nonlinear parameters 0^ for the 
analysis. The additional information needed to set up the above adaptive 
design procedure is hence minimal. Obviously, any other strategy that al- 
lows one to use a set of candidate dose-response models might also be used 
for the final analysis. 

3.4. Technical details. In this section we provide further details of the 
algorithm presented above. 

For the calculation of the one- and two-dimensional integrals in (3.5) we 
used quasi-uniformly distributed point sets based on good lattice points 
ui,...,u„, where Uj G [0,1]'^ and d is the dimension of 6^; see Fang and 
Wang (1994) for details on the construction of such integration grids. Let 
■K{6m\y) =p(y|^m)'^) ><p(^ml"^) denote the integrand from (3.5) and let b^ 
and b„ denote the vector of lower and upper bounds for 0^. One first trans- 
forms the good lattice points to obtain u* = Uj(b„ — hi) + hi for i = 1, . . . , n, 

and then approximates the integral (3.5) by Ylj=i{buj — hj) Yll=i '^{^t ll/)/''^- 
This approach also allows one to calculate the approximate maximum in the 
subsequent optimization step by using the grid point u* corresponding to 
maxj7r(u*|y). We found that using a grid of size 100 in the one-dimensional 
case and the good lattice point set of size 1597 in the two-dimensional case 
[Fang and Wang (1994)] provide sufficiently reliable and computationally 
efficient results (both for integration and optimization). 

The optimization in (3.4) is a constrained optimization problem because 
the weights Wncxt lie in the {k — l)-dimensional probability simplex. A sim- 
ple but efficient approach to perform the optimization is to use a mapping 
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]^fc-i |__j^ gk g^j^^ then to employ a standard unconstrained optimizer, as de- 
scribed in Atkinson, Donev and Tobias (2007), page 131. To account for 
the already allocated patients until an interim analysis, one can optimize 
^( ""'Aftif+M.^r'" ) '^^^^ respect to Wnext € S''- Due to potential multiple op- 
tima in the design surface, one cannot be sure whether indeed an optimal 
design has been found by the optimizer. We thus propose using lower bounds 
of the resulting relative efficiencies based on the underlying geometry of the 
optimization problem. 

To be precise, suppose that the vector w* has been found by the opti- 
mizer. The following result gives a lower bound on r(w*) = q,T°t\ S [0, 1], 
where Wopt is the (unknown) true optimal design at the end of the next 
stage, accounting for the patients allocated until the current interim analy- 
sis. A proof of the result is given in the Appendix. 

Theorem 3.1. A design w with Cm{Om) = '^bm{Om) G Range(Mm(0m, 
w)), m = l,...,M, minimizes ^(w) with respect to Wnext; where w = 
°N _"jv'''^"°'"' ; if and only if there exist generalized inverses Gi,. . . ,Gm 
o/Mj„(0miw); such that the inequality 



J2m=l(^rnCm{dm)G'^Mm{0m,Wr,ext)GmCm{0m)/c^{0)GmCm{d) 



is satisfied for all d € {di, . . . , dk}- Moreover, the efficiency of any design w 
can he hounded from helow hy 

('•^^ ^("^) ^ F(^ ^ 7^' 

where 7 = ^^^^^^^..1 ' ^*(w) = minGi,...,G„ maxdg{rfj^...^dj /i(d, w), 



F(w,7) 



M 

T 



l + (l-7) min max y^ amcLiGm)G 



m=l 



X (Mm(0m,v) -Mm(0m,Wncxt)) 



X GmCm{d)/icl,{9m)GmCm{0m)) 



-1 



and the minimum is taken over all generalized inverses Gi, . . . ,Gm of the 
matrices Mi(0i,w), . . . ,Mm(^m, w). 

When the matrices Mi(0i,w), . . . ,Mjv/(0a/,w) are invertible, /i*(w) is 
just the maximum of h{d,w) over the k doses and straightforward to calcu- 
late [and so is the lower bound on r(w)]. This lower bound is useful in several 
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respects: we do not need to know the actual optimal design Wopt in order to 
calculate the lower bound. If the lower bound for our calculated design w* 
is equal to 1, we know that w* is the optimal design. Otherwise, we have 
a conservative estimate on how much percent off one would be when using w. 
The bound based on A;*(w,7) is sharper, however harder to implement. 

If one does not use a fully Bayesian approach for the final analysis, one 
typically has to fit nonlinear regression models to the data. When there are 
only few doses available, as it is often the case in drug development prac- 
tice, calculating the ML estimate may be difficult. One way to simplify the 
problem is by exploiting the fact that 6q and di enter the model function 
linearly in (3.1). We thus apply the nonlinear optimization only on the non- 
linear parameters 0^, similar in spirit to Golub and Pereyra (2003). Using 
the Frisch-Waugh-Lovell theorem [Baltagi (2008), Chapter 7], we can recal- 
culate the residual sum of squares efficiently, without the need to solve the 
full least squares problems in each iteration of the nonlinear optimization 
(this effect becomes even more important when there are additional linear 
covariates in the model equation, such as gender, baseline values, etc.). In ad- 
dition, we impose bounds on the nonlinear parameters 0„^ to guarantee the 
existence of the least squares estimate [Seber and Wild (1989), Chapter 12]. 
As mentioned in Section 3.2.1, such bounds are not a severe restriction in 
practice and ensure that the optimization problem is well posed. 

4. Asthma study revisited. In this section we revisit the asthma case 
study from Section 2 and address the four open design questions using the 
proposed methodology from Section 3. To this end, we investigated in an 
extensive simulation study the operating characteristics for different design 
options and parameter configurations. 

4.1. Design of simulation study. We generated normally distributed ob- 
servations according to the dose-response models given in Table 1 with 
a = 350 ml. To investigate the robustness of the proposed methods, we also 
simulated from a linear model (with baseline 100 ml and maximum effect 
300 ml) that was not included in the candidate model set. The total sample 
size was fixed at 300 (constraint imposed by the clinical team). To evaluate 
the benefit of including additional doses, we compared two design options, 
one with the four active doses 2.5, 10, 20, 50 (plus placebo) and another one 
with the seven active doses 0.5, 1, 2.5, 5, 10, 20 and 50 (plus placebo). In 
addition, we evaluated the benefit of additional interim analyses by varying 
their number from (= no interim analysis) to 9, where the interim looks 
were chosen equally spaced in time. In all cases, we assumed a balanced first 
stage design. The designs from the second stage onward were determined 
using the observed data according to the algorithm from Section 3.3. When 
the MED estimate did not exist for certain models at an interim analysis. 
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they were removed from the model set for the purpose of design calculation 
and the model probabilities were reweighted accordingly. When the MED 
estimate did not exist for any model, a balanced allocation was used for the 
next cohort of patients. 

For the final analysis we employed the MCP-Mod procedure from Bretz, 
Pinheiro and Branson (2005). A potential dose-response signal was assessed 
using model-based multiple contrast tests based on the candidate model 
set from Table 1. Subsequently, if there were significant models, the dose- 
response model with lowest Akaike Information Criterion (AIC) among the 
significant models was chosen to estimate the MED. 

The methodology from Section 3 was applied with uniform prior proba- 
bilities for the different models. We further assumed a priori distributions 
with mean 100 and variance 100,000 for the placebo effect and mean 300 and 
variance 100,000 for the maximum treatment effect, which were then trans- 
formed into the linear parameters for all dose-response models. The mode 
of the marginal distribution for a"^ was chosen as 350^ with z/ = 4, result- 
ing in an infinite variance. For the nonlinear parameters we assumed beta 
distributions (or products thereof) with mode equal to the values specified 
in Table 1 and S" = 3. The parameter bounds were chosen to ensure that all 
reasonable dose-response shapes remained included within the bounds. That 
is, we chose O2 G [0.05, 75] for the Emax models, (^2, ^3) G [0.5, 4] x [0.5, 4] for 
the beta models and (^2,^3) G [0.05,75] x [0.5,25] for the the logistic model. 
For each scenario we used 5000 simulation runs. 

4.2. Simulation results. For the chosen standard deviation of o" = 350 ml, 
the power of the MCP-Mod procedure to detect a dose-response signal was 
almost always close to 1. Thus, the MCP-Mod procedure was essentially 
reduced to choosing the nonlinear model with lowest AIC value under the 
constraint that only models with significant contrast test statistics were 
included in the model selection step. Simulations with a values larger than 
350 ml indicated that the power quickly dropped to lower levels (results not 
shown here), although the estimation results remained qualitatively similar 
to the ones shown below. 

In Figure 3 we display, for each simulation scenario, the mean absolute 
estimation error for the MED against the number of interim analyses. In 
all scenarios one observes a benefit from adapting, while most of the im- 
provement is already achieved after 1, 2 or 4 interim analyses. The largest 
relative improvement (comparing no adaptations vs 9 adaptations) can be 
observed for the Logistici and the Beta model scenarios, particularly in the 
case of 7 active doses. The worst relative improvement can be observed for 
the Logistic2 scenario, where the overall largest absolute estimation error 
can be observed. This is not surprising, because even when adapting one 
cannot achieve a good design for this model, as there are no doses available 
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Fig. 3. Mean absolute estimation error for MED estimation. 

for administration in the interval (20, 50) containing the MED; see also Fig- 
ure 1. It is remarkable to see that adaptation also works in the linear model 
scenario, although the linear model is not included in the candidate model 
set. It seems that other models in the candidate set are able to capture the 
shape of a linear model reasonably well. 

The comparison between 4 and 7 active doses is not entirely clear. If no 
interim analyses are performed, it seems that the design with a balanced 
allocation across the 4 active doses is slightly better than the design with 
a balanced allocation across all 7 active doses. If one decides to adapt, how- 
ever, it seems beneficial in some cases to have more doses available, partic- 
ularly if many interim analyses are performed, while in other cases 4 active 
doses are sufficient. 

To illustrate how adaptation changes the allocation of patients to the 
different doses, we display in Figure 4 the average patient allocations for 
the Emax2 model after 1, 2, 4 and 9 interim analyses and with 7 available 
active doses. The adaptive design tends to allocate more patients both on 
placebo and nearby the actual MED. This is intuitively plausible, as the 
MED estimate depends on the precision of the estimated placebo effect 
as well as of the estimated function /(•) around the true MED. It also 
follows from Figure 4 that for a large number of interim analyses the overall 
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model. Last panel: locally MED-optimal design for the Emax2 model with true MED = 7. 7. 



allocation is close to the one under a locally optimal design for the Eniax2 
model, with the variability in the allocations due to the uncertainty both 
in estimating the correct model and the model parameters at the interim 
analysis. Similar conclusions also hold for other models than the Emax2 
model (not reported here). 

We now investigate to which extent the precision gain observed in Figure 3 
translates into sample size savings when performing an adaptive design. In 
other words, how many additional patients are required for a nonadaptive, 
balanced design to achieve a similar estimation error as with an adaptive 
design using 300 patients. We again considered the Emax2 model and it- 
erated the total sample size until the mean absolute estimation error was 
approximately 4 (which is the mean absolute estimation error obtained af- 
ter 9 interim analyses, as seen in Figure 3). For both design options with 4 
and 7 active doses, this was achieved after roughly 500 patients. Thus, using 
a nonadaptive, balanced design, one would need 200 additional patients to 
achieve a similar precision in MED estimation as compared to an adaptive 
design using 300 patients. 

The adaptive design benefits observed so far depend on several input pa- 
rameters, such as the starting design for the first stage. One may argue 
that starting with a bad design that allocates patients at the "wrong" doses 
may be improved by adapting at one or more interim looks. On the other 
hand, starting with a good design may lead to adaptations following ran- 
dom noise at the interim analyses. To illustrate this effect, we report the 
results for the simulations under the Logistici model (similar results were 
also obtained for other models and scenarios, but are not reported here). We 
used four different starting designs. We used w = (0.35,0.03,0.22,0.35,0.05) 
and w = (0.35,0.02,0.02,0.02,0.02,0.20,0.30,0.07) as good starting designs 
with 4 and 7 active doses, respectively. These designs work well because they 
allocate patients on placebo and around the MED, while keeping some mass 
on the remaining doses. In addition, we used w = (0.1,0.3,0.05,0.05,0.5) 
and w = (0.1,0.2,0.22,0.02,0.02,0.02,0.02,0.4)' as bad starting designs, as 
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Fig. 5. Mean absolute estimation error for MED estimation, under the the Logistic^ 
model for good and had starting designs. 

they have relatively few patients on placebo and around the MED. It follows 
from Figure 5 that substantial improvements are possible when using bad 
starting designs. On the other hand, for good starting designs no benefit is 
achieved by adapting and the performance may even deteriorate, because 
the possibility of adapting may lead one to deviate from the already good 
starting design. In practice, one does not know whether an employed design 
is good or bad, but one should keep in mind the possibility that adaptive 
designs will not always improve upon the initial design. 

To further investigate the robustness of the proposed methods, we re- 
peated the simulation study from Figure 3 by increasing the standard de- 
viation to 450 ml and 700 ml. The overall results remain similar, but with 
increased absolute estimation errors. However, the relative benefit of 9 in- 
terim analyses vs. no adaptation decreases slightly. Due to the larger noise, 
one obtains less reliable information at an interim analysis and one may end 
up with a worse design for the next stage. We also investigated the effect of 
prior misspecification. For this purpose we misspecified the prior means or 
prior modes by adding or subtracting 20% of the true value, but leaving the 
variability (variance of baseline and maximum effect and the value S for the 
beta distribution) unchanged as in the original simulations. The results are 
largely identical to those presented in Figure 3, indicating that the proposed 
methods are robust under moderate prior misspecifications. 



4.3. Conclusions for asthma study. Many more simulations than pre- 
sented above were conducted at the planning stage of the asthma study to 
address the four questions stated in Section 2. Regarding question (A), it 
was felt that the potential benefits of conducting an adaptive design (more 
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precise MED estimation) outweighed the additional logistical requirements, 
especially in view of the perceived sample size gain of 100-200 patients when 
compared to a fixed-sample study designed to achieve a similar precision. For 
question (B) it was decided to have one interim analysis: based on Figure 3 
and other simulation results, the potential further reduction of the mean 
absolute estimation error with two or more interim analyses was perceived 
as too small to justify the additional logistical complexity. 

For similar reasons, it was decided against having all seven actives doses 
from the beginning on question (C). Instead, 150 patients ought to be allo- 
cated equally across the four active doses 2.5, 10, 20, 50 (plus placebo) in 
the first stage. Once the interim results are available and analyzed with the 
methods from Section 3, however, patients could be allocated to all seven 
active doses (or a subset thereof) in the second stage. For practical rea- 
sons, the clinical team decided to incorporate constraints on the minimum 
number of patients allocated per dose in the second stage: if the algorithm 
would allocate less than 5% of the patients on a certain dose, that dose 
would be dropped altogether and the corresponding patients reallocated to 
the remaining doses. 

5. Discussion. Motivated by a dose finding study in moderate persistent 
asthma, we described a response-adaptive approach that addresses common 
challenges encountered in dose-finding studies: multiple objectives, model 
uncertainty, and large variability. When planning an adaptive dose-finding 
design it is important to realize that it may not always be better than 
a nonadaptive design. It is necessary to employ a factored view, as many 
parameters may impact the performance of a study design. Often, an unbal- 
anced fixed-sample design derived from optimal design theory might already 
provide benefit over a balanced fixed-sample design and adaptation may not 
bring further advantages, particularly if the variability is large (which is 
common in practice). Thus, adaptive designs are promising in situations 
where the initial design is not good and/or interim parameter estimates 
have low variability. In practice, one never knows how good the initial de- 
sign will be, before trial start, and adaptive designs may guard against bad 
initial designs. However, the benefits of adaptive dose-finding designs have 
to be balanced against the increased logistical requirements to implement 
processes for repeated data collection, cleaning and analyses, to maintain 
trial integrity and validity, and to overcome potential challenges in drug 
manufacturing and supply. 

In this paper we focused on designs based on the compound optimal- 
ity criterion (3.4) to address model uncertainty and to minimize the vari- 
ance of MED. The criterion depends on the parameters of the different 
dose-response models as well as on the model probabilities and we used 
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a Bayesian approach to continuously update parameter values and model 
probabilities based on the information accrued in the trial. The approach 
was implemented based on optimization and numerical quadrature, so that 
computationally intensive Markov chain Monte Carlo techniques could be 
avoided. Computational efficiency is of extreme importance, as the frequen- 
tist operating characteristics of any adaptive design methodology needs to 
be evaluated in extensive simulations under multiple scenarios. 

The proposed method can be extended immediately if alternative opti- 
mality criteria [such as EDp- or D-optimal designs, see Dette et al. (2010)] 
or mixtures thereof are of interest. Alternatively, optimal discrimination de- 
signs could be applied that allow one to differentiate among several candidate 
nonlinear regression models [Atkinson and Fedorov (1975), Dette and Titoff 
(2009)]. It would be interesting to address multiple objectives by considering 
different optimality criteria at different stages, such as using a model dis- 
crimination design in earlier stages, and MED-optimal design in later stages. 
This will be investigated in future research, but see Dragalin et al. (2010) 
for initial results. 

The R functions used for the simulations are available with the 
DoseFinding R package [see Bornkamp, Pinheiro and Bretz (2010)]. 

APPENDIX: PROOF OF THEOREM 3.1 

Obviously the first part of the theorem follows from the lower bound (3.7) 
on the efficiency. For a proof of (3.7) let 7 = A''oid/(A''oid + A^ncxt) G (0, 1) and 
note that the total information of the experiment in the mth model is given 
by 

(A.l) Mm(0, Wold, Wncxt) = 7Mm(0m, Wold) + (1 - 7)Mm(0m, Wnext), 

where we collect in the vector 6 = (61, . . . , Om) the parameters of the differ- 
ent models. Define a block diagonal matrix by 

M(0m,Wold,Wncxt) 
(A.2) 

= diag(Mi(0.m, Wold, Wncxt), • • • ,Mm(0m, Wold, Wnext)) 

(all other entries in this matrix are 0) and, similarly, 

K = diag(ci(0i), . . . , cm{Om)), 

where the vector CmiO) is given by CmiO) = '^hm{0),m = 1, . . . , M. For a de- 
sign Wnext, such that Cm{9) G Range (Mm (0m, Wold, Wnext)) (m = 1, . . . ,M), 
we consider the information matrix 

Ci^(M(6',Wold,Wncxt)) = (i^^M"(6>,Wold,Wnext)i^)"^ 

= diag((cf (0i)Mj"(0i, Wold, Wnext)ci(0l))"\ . . . , 

(cM(0M)M^(0A/,Wold,Wncxt)cA/(0A/))~"^)- 



20 BORNKAMP, BRETZ, DETTE AND PINHEIRO 

Note that the optimal design maximizes 

^-^Wnext) = ^°'^+^"^-^ . cI>,(Cx(M(0,WoM, W„ext))) 

where the last identity defines the criterion <1>q, and we have used the notation 
<^a(diag(Ai, . . . , Am)) =11™=!-^™"- Now according to Theorem 1 in Dette 
(1996), a lower bound for the efficiency of the design Wncxt 

V'-Hw) $„(Cfc(M(0,Wold,Wnext))) 



r(w) 



-0 ^Wopt) max^ggfc$<^(Ci^(M(0,Woid,v))) 
is obtained as 

minmaxtr{GE:Ci^(M(0, Wqm, Wnoxt)) 
. G AeA 

X ECi^(M(0,Wold,Wnext))i^^G^A}]"\ 

where the minimum is taken over the set of all generalized inverses of the 
matrix M(0, Woid,Wnext) and the set A is defined by 

^={M(0,Woid,v)|vGS^-} 

and the matrix E is given by 

E = diag(aicf(0)M^(0i,v^oid,v)ci(0),...,aMCM(0)Mj^"_^(0M,Woid,v)cm)- 
Therefore, observing the identity 

Mm{Om, Wold, v) = Mm{9m, Wold, Wncxt) + (1 -7)(M(0m, v) - M(0, V^next)), 

we obtain 

M 

1 + (1 - 7) minmax V] amC^(0m)G^(Mm(0m, v) - M.m(0.m, Wncxt)) 
Om lies'-' , 

m=l 

-1 

X GmCm{0)/ic^{e„r)G„rCmi9m)) 



> 



mm max > a, ~ "* 



idmi'^^ 9rn)GmCm{6m)) 



Gu...,G,nde{di,...,dk}^-^ '" cJn{Orn)GrnCrn{Om) 



E Cm{9m)G^Mm{6m,'^ncxt)GmCm{9m 
"'" c'^iO )G C (6 ) 



RESPONSE-ADAPTIVE DOSE-FINDING UNDER MODEL UNCERTAINTY 21 
where we have used the inequahty 

[i+(i-7)(^-i?r^> 

for A> B >0, (1 — j)B <1 and standard arguments in design theory. 
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